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1 We present an analysis of the magnetic response of a mesoscopic superconductor, i.e. a system 

of sizes comparable to the coherence length and to the London penetration depth. Our approach 
is based on special properties of the two dimensional Ginzburg-Landau equations, satisfied at the 
£Nj \ dual point (ft = -j^)- Closed expressions for the free energy and the magnetization of the supercon- 

ductor are derived. A perturbative analysis in the vicinity of the dual point allows us to take into 
account vortex interactions, using a new scaling result for the free energy. In order to characterize 
the vortex/current interactions, we study vortex configurations that are out of thermodynamical 
equilibrium. Our predictions agree with the results of recent experiments performed on mesoscopic 
aluminium disks. 
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Cil' I. INTRODUCTION 

The ability to detect and manipulate vortices with great sensivity in systems of small size such as mesoscopic 
superconductors [Q or atomic condensates Q has generated an outgrow of interest in the mechanism of creation and 
annihilation of vortices and in the study of stable and metastable vortex configurations. In particular, recent advances 
I | in the technique of Hall magnetometry 101 have allowed to measure the magnetization of small superconducting 
samples containing only a few vortices JljQ. These experiments are conducted on aluminium disks well below the 
superconducting transition temperature, whereas previous measurements were performed only in the vicinity of the 
normal/superconductor phase boundary (^||. Besides, the magnetization measurements in Jjj] are carried out on an 
individual disk and not on an ensemble of disks as in |6| . The radius R and the thickness d of the sample used in 
the experiments are comparable to the superconducting characteristic lengths, i.e. the London penetration length 
(A = 70nm) and the coherence length (£ =250nm). Such a sample can neither be considered to be macroscopic, nor 
microscopic. The system falls, rather, in a mesoscopic regime where surface effects are of the same order of magnitude 
as the bulk effects. Thus, the magnetic response of a mesoscopic superconducting disk to an applied field depends 
| strongly on its size and is very different from that of a macroscopic superconductor. When the radius R of the sample 
I/"") . is much smaller than the coherence length £, no vortex can nucleate, the normal/superconductor phase transition is 
second order and the magnetization M, as a function of the external applied field H e , has a non-linear behaviour (non- 
linear Meissner effect If R is comparable to £, the superconducting phase transition is first order and a bistable 
hysteresis region appears in the M — H e curve. For R greater than £, the phase transition is again second order, and 
when the applied field exceeds a critical value H±, the magnetization curve exhibits a series of discontinuous jumps 
corresponding to the successive entry of vortices into the sample. This qualitative interpretation is supported, at least 
for low applied magnetic fields, by the periodicity of the jumps which corresponds to the entrance of an additional 
superconducting quantum of flux into the disk. For larger fields, or equivalently for higher density of vortices, both 
the period and the height of the jumps become smaller, a behaviour related to the interactions between the vortices 
O , and to transitions between stable vortex configurations, with the same number of vortices. 

O ■ The magnetization shows also a hysteretic behaviour depending on the direction of the field sweep, due to the 
, presence of a confining energy barrier (the absence of remanent magnetization precludes pinning effects). In some 
metastable states, the sample may exhibit even a paramagnetic response whereas in thermodynamic equilibrium 
a superconductor is diamagnetic. 

These experimental results have led to a renewed interest in the theory of mesoscopic superconductors. Numerical 
computations have shown that the phenomenological Ginzburg-Landau theory is well suited to describe a supercon- 
ducting sample in the mesoscopic regime, even far from the critical temperature. These works have revealed physical 
phenomena that play an important role in such systems (for a review see Q), like: the role of surface barriers for vor- 
tex nucleation and hysteresis |9p|i~l|]; the interplay between vortex/vortex and vortex/edge interactions that explains 
vortex structures in mesoscopic disks 10 .[J; the transition between a giant multiple vortex state and a state with 



several vortices carrying a unit quantum of flux [ |12| 

The Ginzburg-Landau free energy of a superconductor involves two fields, the (complex) order parameter tp = \tp\e lx 
and the vector potential A. The minimization of this free energy leads to a set of two coupled non-linear partial 
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differential equations for ip and A, involving the two characteristic lengths A and £. But the solutions depend only on 
one relevant number, the phenomenological Ginzburg-Landau parameter k defined by 

«=j. (1) 

A macroscopic superconductor is said to be of type I if n < and of type II if k > A- . A macroscopic superconductor 
of Type II admits a stable Abrikosov vortex lattice phase when the applied field H e lies between the first penetration 
field and the upper critical field [fT3"| . For aluminium, k is smaller than -^=, hence a macroscopic sample of Al is a 
type I superconductor. 

Analytical studies of the Ginzburg-Landau equations in two-dimensional systems require the use of various approx- 
imations since, in general, exact solutions can not be found due to the non-linearity. One approach is to linearize the 
equations assuming \^)\ <C 1, and to decouple them by supposing that the magnetic field B in the sample is equal 
to the applied field H e . This approach describes correctly the superconducting/normal phase boundary ]^||[l4j-[l6[ , 
but fails to explain the behaviour of the sample deep inside the superconducting state. For example, in the linearized 
theory, all the vortices are at the center of the disk |l6| and therefore one can not study the role of surface barriers, 
the interaction between vortices, and the fragmentation of a giant vortex into unit vortices. Besides, the critical fields 
corresponding to the successive entrance of vortices into the sample do not scale correctly with the size of the system 
(e.g. experimentally, the entrance field H\ of the first vortex scales as whereas the linear theory predicts a Rr 2 
dependence). Of course, in the vicinity of the upper critical field ]T^ ] the linearized theory agrees quantitatively with 
the experimental results. 

A second approach is to use the London equation which can be derived from the Ginzburg-Landau equations by 
supposing that |^| = 1 everywhere except on a finite number of isolated points, called vortices, where \tj)\ = 0. 
London's equation is valid rigorously when the parameter k goes to infinity, i.e. for extreme type II superconductors 
in which vortices are indeed point-like. Many theoretical results have been derived from the London equation, such 
as discrete nucleation of flux lines in a thin cylinder or in a thin disk p^ , p0[ |, the existence of surface energy 

barriers [ p^|j2l| , and the computation of polygonal ring configurations of vortices in finite samples |22] , p3]] . However, 
when k — > oo, the minimum energy is obtained for one flux quantum per vortex [^4|,^5| and vortices have a hard-core 
repulsive interaction impeding the formation of a giant vortex state. Moreover, the surface energy barriers calculated 
from London's equation are quantitatively different from those obtained by numerically solving the Ginzburg-Landau 
equations [p6| . In fact, the experimental conditions are far off the London limit, although thin Al disks are likely to 
have an effective k greater than its measured value Q of 0.28 (in a thin disk, one can argue, following [p^| , that the 
effective London length is of the order of A 2 /d, and this results in a higher value of k). 

We follow another approach, less explored in the literature, based on an exact result for the two-dimensional 
Ginzburg-Landau equations. In an infinite plane reduce to first order differential equations that can be decoupled 
when the parameter k takes the special value 1/V2, called the dual point 24 2^J|^]. At that point, the free energy 



is a topological invariant of the system. In ]28[ , we generalized this method to a finite domain with boundaries; this 
enabled us to classify solutions with different number of vortices and to derive analytical expressions for the free 
energy and the magnetization of a mesoscopic disk as a function of the applied field. Our results agreed qualitatively 
with the experimental data, and even quantitatively when the number of vortices in the system is low. However, some 
important features such as the non-linear Meissner effect in a fractional fluxoid disk, the variation of the amplitude 
and the period of the jumps in the M — H e curve could not be described. Moreover, in plj , we discussed only the 
case where R is much larger than £ and did not obtain the different regimes of the magnetization curve when the 
ratio is varied. 

In this paper, we study the Ginzburg-Landau free energy T not only at the dual point k = l/v2 but also in its 
vicinity where vortices start to interact weakly p9j . Taking into account non- linear effects, our calculations describe 
the magnetic response of the sample as its size changes, providing an understanding of the non-linear Meissner effect 
and of the multi-vortex state. We shall also study non-equilibrium vortex configuration in order to determine the 
interaction between a vortex and edge currents. 

The plan of this paper goes as follows: in section 2, some basic features of the Ginzburg-Landau theory of super- 
conductivity are recalled. In section 3, after studying the case of an infinite system, we generalize the Bogomol'nyi's 
approach to a finite size superconductor and calculate its free energy at the dual point. This result is applied to 
an infinite cylinder in section 4. The case of a mesoscopic disk is studied in section 5 and magnetization curves are 
obtained for systems of different sizes. In section 6, we obtain the free energy and the magnetization of a cylindrically 
symmetric system when k is close to the dual point. The surface energy barrier for a one vortex state out of ther- 
modynamic equilibrium is calculated in section 7. In the last section we discuss our results and suggest some further 
generalizations. Some mathematical details are included in the two appendices. 
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II. THE GINZBURG-LANDAU THEORY OF SUPERCONDUCTIVITY 



We recall here some basic features of the Ginzburg-Landau theory and define our notations. The order parameter 
ip = \ip\e lx is a complex number and the potential vector A satisfies V x A = B, where B is the local magnetic 
induction. The two characteristic lengths A and £ appear as phenomenological parameters. In this work, we measure 
lengths in units of A-\/2, the magnetic field in units of 4 ^° 2 and the vector potential in units of A where the flux 

quantum <f>o is given by (f>o = |f . The Ginzburg-Landau free energy defined as the difference of the free energies 

T = Fs(B) — Fs(0), is measured in units of ^ where H c the thermodynamic field satisfies H c = \/2k J'% . In these 
units, T is given by 

T = /J |S|2 + ^ ^ + l( ^ ~ ' (2) 

where the integration is over the superconducting domain Q. The Ginzburg-Landau equations that minimize 
become 

- (V - iAfi) = 2k 2 V(1 ~ W) (3) 
V x B = 2j (4) 

Equation (||) is the Maxwell- Ampere equation with a current density j = Im(ip*Vip) — \ip\ 2 A related to the superfluid 
velocity v s by 



J_ 



UJ2=VX-A (5) 



Outside the superconducting sample, ip = 0. The boundary condition on the surface of the superconductor is obtained 
by requiring that the normal component of the current density vanishes (superconductor/insulator boundary condition 
II): 

(V - - (6) 

here n is the unit vector normal at each point to the surface of the superconductor. 

J 

function the circulation of the London fluxoid along a closed contour C is quantized ]13|,[30| 



The London fluxoid is the quantity ( | ■ + A] , that is identical to Vx- Since x is t ne phase of the univalued 
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c {r^+A).dl = jVx-dl = 2nn (7) 

The integer n is the winding number of the phase of the system along the contour C and is a topological characteristic 
of the system. 

In this study, the superconducting sample is either an infinite cylinder or a thin disk, with cross-section of radius R, 
placed in an external magnetic field parallel to its axis. Since R is an important parameter, we define the dimensionless 
quantity: 

«=— (8) 

a is supposed to be small compared to 1 (typically a ~ 1/10 in the experiments) unless stated otherwise. The flux 
created by the external and uniform magnetic field H e (expressed in units of j^p-) through the cross section ttR 2 of 
the sample is equal to irR 2 H e = y^'Po- The flux <f> e , in units of the flux quantum <fio, is thus given by: 



: 4ttA2 2a? 

fie = 



(9) 



We emphasize that, in the units we have chosen, the flux <f>b of a magnetic field B through a surface f2 is obtained via 
the following formula: 
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<t>b = ^- I B.dS =^-l A.dl (10) 
^ Jn 2tt Jan 

An extra factor 1/2% appears here because B is given in units of 2 , the surface in units of 2A 2 and the flux in units 

of (j) . 

Since we are studying a superconductor in an applied external field, the relevant thermodynamic potential is the 
Gibbs free energy G obtained from F via a Legendre transformation: 

G = F -H e B = F — H e 2n<j> b = F — Aira 2 <f> e (/> b (11) 
Jn 

In a normal sample, tp = and B — H e . Therefore, the Gibbs free energy Gjv of a normal sample is given by: 

Gn = F N - H e I B — Fn — 2na 2 <p 2 e (12) 
Jn 

At thermodynamic equilibrium, the superconductor selects the state of minimal Gibbs free energy. The quantity that 
we are interested in, and which is measured in experiments, is the magnetization M of the superconductor due to the 
applied field given by 4irM = B — H e . It is obtained, at thermodynamic equilibrium and up to a constant equal to 
the superconducting condensation energy, from the difference of the (dimensionless) Gibbs energies 

Q = G S -G N = T + 2na 2 (g - Ana 2 ^ (13) 

using the thermodynamic relation |L3] : 



M =2 L § ^ 

2lT OQ e 



III. FREE ENERGY OF A SUPERCONDUCTOR AT THE DUAL POINT 

We now study the particular case of the dual point, defined by k = For this value of the Ginzburg-Landau 
parameter, the free energy (0) of a two dimensional domain O can be written as p5|J2S]: 



T = \ (\{B-l + H 2 ) 2 + \V^A + I (j + A).dl (15) 
Jn \ z J Jan 

where the operator T> is defined as T> — d x + id y — i(A x + iA y ) and the second integral is over the boundary of the 
domain Q. 

A. The case of an infinite system 

If we suppose that the domain Q is infinite and superconducting at large distances pj3], i.e. — > 1 at infinity, 
then the boundary integral in (15) is identical to the fluxoid. Using the quantization property (0), we obtain 

? = 27TB + J Q(B - 1 + \yj\ 2 ) 2 + (16) 

The free energy is thus minimum when Bogomol'nyi equations p5[ are satisfied, that is when, 

Vtp = (17) 
B = l -M 2 (18) 

Thus, the total free energy results only from the boundary term in (|l^) and is a purely topological number: 

T = 2im . (19) 

The free energy is proportional to the number of vortices: at the dual point, vortices do not interact with each other 

ppi. 
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B. Finite size systems 



In a finite system with boundaries, vortices do not interact with each other at the dual point but they are repelled 
by the edge currents. Therefore, at thermodynamic equilibrium, all vortices collapse into a giant vortex state. Since 
the superconductor under discussion has a circular cross-section, this giant vortex (or multi-vortex) is located at the 
center and the system is invariant under cylindrical symmetry. In a finite size mesoscopic superconductor at the dual 



point, the boundary integral, in (15), can not be identified with the fluxoid because \ip\ is in general different from 1 
on the boundary of the system. This quantity is no more a topological integer but a continuously varyin g re al number. 
The two terms of ( |l5j ) can not, therefore, be minimized separately to obtain the optimal free energy. In p8| , we found 
a method to circumvent this difficulty: if the system is invariant under cylindrical symmetry, i.e. all the vortices are 
at the center of the disk, then the current density has only an azimuthal component jg. The current jg has opposite 
signs near the center (where the vortex is located) and at the edge of the disk (where Meissner currents oppose the 
penetration of the external field). Hence, there exists a circle F on which jg vanishes ^8|. Along T, we have 

j + A = T^j2 + A = Vx and therefore j> (j + A).dl = 2nn (20) 

The domain f2 can thus be divided into two sub-regions, f7 = f2i U such that the boundary between f2i and Q2 
is the circle T. By convention, we call Oi the bulk and the annular ring Q2 the boundary region (see figure 1). 

Numerical solutions of the Ginzburg-Landau equations in a two dimensional superconductor, with cylindrical sym- 
metry, clearly show the separation of the sample cross-section into two distinct subdomains. In figure 2, we have 
plotted the order parameter and the magnetic field in a cylinder of radius R = 10A\/2 with one vortex at the center. 
These two quantities vary only near the center and near the edge: there is a whole intermediate region in which l^l 
and B remain almost constant. When the system is large enough, these constant values are, to an excellent precision, 
identical to the asymptotic values of \tp\ and B in an infinite system. The current vanishes for a value of r for which 
= 0, and this determines the radius of the circle T. In figure 2, this corresponds approximately to r ~ 5.5, though 
practically, the circle T can be placed anywhere in the saturation region where the current is infinitesimally small. 

Although the Bogomol'nyi equations ( |l7| , fl8l ) do not minimize the Ginzburg-Landau free energ y [jjll , we notice that 
the behaviour of \4>\ and B in the bulk subdomain Q,\ is still given by the relation B = 1 — |-0| 2 ,(|l8|) which, as shown 
in figure 3, represents indeed an excellent approximation. Thus, using ( |l5| ) and (^0|), we conclude that at the dual 
point the free energy of Oi can be calculated as that of an infinite domain, namely 

F(Qi) = 2nn (21) 

We also emphasize that the flux in f2i is quantized and one has: 



1 

2^ 



B = n (22) 

Hi 



To calculate T(£l?), the identity ( |15| ) valid at the dual point is of no use anymore, since \ip\ is in general different 
from 1 at the boundary, and the boundary integral in ( |l5| ) can not be identified to the fluxoid. We therefore have to 
go back to the definition @ of the Ginzburg-Landau free energy which becomes at the dual point: 

*W = J n % + (VM) 2 + M 2 |V X - A\ 2 + (1 ~!f |2)2 (23) 

The assumption of cylindrical symmetry implies that x — n ® where is the polar angle and n the number of vortices 
present at the center of the disk. Examining again figure 2, we observe that in O2 , the order parameter and the 
magnetic field vary from their values on the edge to their saturation values over a region of width 5, which is of order 



1 in units of Av2[| The length S therefore represents the typical distance over which the integrand in (23) has a non 
negligible value. 

With the help of this observation, we shall estimate J- (^2) using a variational Ansatz: we shall consider that the 
modulus of the order parameter has a constant value ipa over a ring of width 5, included in SI2 and that A and B decay 



1 Indeed, one has for a thick system S ~ A at the dual point. For a thin film of thickness d, S ~ A 2 jd in the London limit [|l9[ 
Since we are considering a mesoscopic regime in which d ~ A, both expressions indicate that <5 is of order 1. 
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exponentially with a characteristic length S from their boundary value to their bulk value. Clearly, our approximation 
will be valid only if the width of CI2 is large enough compared to 1. We first remark that our Ansatz is compatible 
with the boundary condition (||), which reduces here to ^ = and that it allows us to neglect the curvature term 
(V|-0|) 2 in (f23|). To evaluate the term proportional to the superfluid velocity v s (r) (||), we first notice that, due to the 
Meissner effect, it decreases from the boundary at r = R with a behaviour well described by 

v s {r)=v s (R)e-^' 5 (24) 

with v s (R) = a(n — 4> b ). To obtain the last equality we used that the boundary value of the vector potential is 
A(R) = a(f>bUo, where <p b is the total flux through the system. Hence, for a constant amplitude iJjq of the order 
parameter, we have 



The magnetic contribution in ( p3| ) is obtained from the typical magnitude B of the magnetic field in Q2 determined 
using (|l0|) and (§|) as 

cf> b = 1 L[B=±[ B + ^( B = n+-B (26) 
2vr J n 2tt J Ql 2tt J n2 a 

Thus, using the fact that B 2 decrease exponentially with a characteristic length 5/2, we estimate the contribution of 
the magnetic energy to ^(^2) as being: 

1 [ B 2 SB 2 a , , , 2 

After substituting this expression into ( p5|) we minimize T(£l-i) with respect to tpo- The optimal variational value of 
■00 is given by: 

2 _ f f-i Vs 2 (i?) = l-^(n-0 b ) 2 if |a(n-0 6 )|<V2 (2R] 
^°~\0 if \a(n- ( f> b )\>V2 { ' 



Inserting these expressions in (Eq), we obtain the variational free energy T{VL2): 



1 ( Av 2 (R) - Bvt(R) if \a{n-<p b )\<j2 

^ { ] ~\Ta + £svl{R) if \a{n-4> b )\>V2 ^ ] 



$ ( 1 
A= — 1 



with A and B defined by 

A = 

2a V 2S 2 

The total free energy of the mesoscopic superconductor containing n vortices, at the dual point, is thus: 

2ir ' \ 5/2a + v 2 (R)/4a6 if \a(n-4> b )\ > V2 v ' 

This energy is the sum of two contributions: 

(i) a bulk term proportional to n which is a topological quantity at the dual point. 

(ii) A boundary term, reminiscent of the well-known 'Little and Parks' free energy fl3|| (this boundary term can be 
given a geometric interpretation in terms of a geodesic curvature p8 32 1). 
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IV. FREE ENERGY AND MAGNETIZATION OF A CYLINDER AT THE DUAL POINT 



We now apply the relations ( |3l"| ) to the simple case of an infinitely long superconducting cylinder of radius R > A, 
lying in an external field H e directed along its axis. There are two contributions to the total flux (pf,- the flux of n 
vortices present at the center of the sample and a fraction of the applied flux </> e localized near the boundary and 
proportional to X/R (due to the Meissner effect). Hence, 

4> b = n + 2acf> e with <f> e = (32) 

The exact numerical coefficient in front of the term a<j) e does not affect the result of our calculation; we take it equal 
to 2, the value obtained in the London limit |jl7|| . The total free energy, using ( |3l| ) and the fact that v s (R) = —2a 2 <fi e , 
is given by 

1 J 4a 4 {A4>1 - 4o 4 B#) if a 2 |0 e | < l/y/2 . . 

^ (n '^ )=n+ U + ^ * i^i>W2 (33) 

Using ( |lT| ) and (|33|), the Gibbs free energy, Q{n, (f) e ), of a cylinder containing n vortices at the dual point is given by: 

^g(n,<t> e ) = n(f-2a 2 e ) + P(4> e ) (34) 



where P(0 e ) is a polynomial in (ft e that does not depend on n. Hence, all the curves G(n, 4* e ) meet at 



(35) 



For values of </> e less than this critical value, the free energy is minimized if there are no vortices. At 4> e = all 
vortices are nucleated simultaneously and the sample becomes normal. This value corresponds to a critical applied 
field H e which is equal to 1 in our units, or restoring the units back, and recalling that k = I/a/2 

^ = ^ = ~^~ (36) 

This is precisely the formula for the thermodynamic critical field of a superconductor |fl3f (which, for a cylindrical 
superconductor with k < , is the same as the upper critical field) . The magnetization M of the cylinder satisfies 
the linear Meissner effect: 

The macroscopic result |fl| is — M = H e ; the finite-size correction to the susceptibility is proportional to i2 _1 . 

Thus, the well-known results for an infinite superconducting cylinder can easily be retrieved from the dual point 
approach. We now proceed to the study of the magnetic response of a thin disk. 



V. A MESOSCOPIC DISK AT THE DUAL POINT 

To modelize the experimental sample of jf|,Q, we consider a mesoscopic disk of thickness d smaller than £ and A. 
Because the disk is very thin, we take the order parameter and the magnetic field to be constant across the thickness 
d of the sample |28|]. This enables us to study the disk as an effective two-dimensional system. However, unlike the 
case of a long cylindrical sample, strong demagnetization effects are present in a thin disk. The value of B near the 
edge of the disk is larger than the applied field H e because geometric demagnetization effects induce a distortion of 
the flux lines |J. Hence the continuity condition B(R) — H e (|32| ) valid for a long cylinder does not apply to describe 
a thin disk. 

In order to find a more suitable choice for the boundary condition for a thin disk, we notice that the higher value 
of the magnetic field at the boundary, a feature which has been obtained from numerical computations [ p3[ , results 
from a demagnetization factor AT close to one, such that jl3| H = in the Meissner phase. The flux lines are 

distorted by the sample and they pile up near the edge of the disk. To describe this, we shall thus take as boundary 
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condition for a thin disk, the expression proposed in |2§[] , which consists in taking the potential- vector at the edge of 
the disk equal to its applied value, i.e. 

A(R) = (j) e au e (38) 

or 

<t>b = 4>e (39) 

Again, this relation does not mean that the field B is uniform and equal to its external strength. A more refined 
value for the boundary condition could have been obtained by using the expression M ~ 1 — ^ j| in the limit d <C R 
of a flat disk. Then, H ~ ^jH e or equivalently 0& ~ -j0 e - But, since (5 ~ d, we shall use for convenience the simpler 
boundary condition given above. 

Substituting v s (R) = a(n — <p e ) in (|3l|), the free energy T(n, <f> e ) of a thin disk containing n vortices is found to be 

Aa 2 (n-<p e ) 2 -Ba^n-^f if a\(n - cj> e )\ < ^2 

± + f s {n-<j> e f if a\(n-cf )e )\>V2 [ ^ } 

and the corresponding Gibbs free energy is obtained using (|l3| ) . In our previous work ps}| , we obtained an expression 
which can be retrieved from ( flO| ) by taking 5=1 and by neglecting the magnetic energy as well as the a 3 term. 
Despite these crude approximations, our analytical results agreed satisfactorily with experimental data, though they 
could not describe neither the behaviour of a disk with a radius smaller than A and £, nor its behaviour when R is 
increased. We apply our present approach to a thin disk with a radius R much smaller than £, and then we consider 
the case R > £. 



2- 



T(n, 4> b ) = n + 



A. Fractional fluxoid disk and Non-Linear Meissner Effect 



We now consider a disk small enough so that no vortices can nucleate i. e. its radius R is less than £ (such a system 
is sometimes called a fractional fluxoid disk [ffl]). If there are no vortices, the domain f^i is empty and Q — f^- Since 
the radius of fl is small with respect to both A and £, we can no more use the expression (|0|) for the free energy, 
but we can assume that the amplitude \ip\ of the order parameter has a uniform value tpo all over the disk and that 
the magnetic field equals the external applied field B — H e . Moreover, in the absence of vortices Vx = and we can 
choose the Landau gauge A(r) = rB/2. Starting from (p3|), and after minimizing the free energy with respect to ipo, 
we find the difference between the free energies of the superconducting and the normal states to be: 



I4 >-T*'i tf 



Q_ 

2tt 

Q 

— = otherwise (41) 
2ir 



From @ we deduce the magnetization M of the sample: 



M=—^- = -(6 e -—4>l) if cub e <V2 



2tt d(j) e 2 

M = otherwise (42) 

The curve representing this magnetization is a cubic. The upper critical field is 4> e — 1/a, i.e. H e oc this scaling 
agrees with the linear analysis of |l6| ] in the limit R £. The transition between the superconducting phase and 
the normal phase is of second order. In figure 4, we plot the relation ji^ ) for —M as a function of the external flux 
(j) e . The dots represent the experimental points obtained from The analytical curve has been scaled so that the 
maximum value of the magnetization and the critical flux coincide with the corresponding experimental data. 



B. Mesoscopic disk with vortices 

We now consider a disk with R > £. The Gibbs free energy difference G(n>, <j> e ) of the disk with n vortices is given 
by (|l3"|). The entrance field H n of the n-th vortex is obtained by solving the equation Q{n,4>e] = Q{n — 1, <fi e ) which, 
using (f40|), reduces to 
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§ = a(l + ^) ((n - 1 - e ) 2 - (n - 0e) 2 ) - ^ ((n - 1 - e ) 4 - (n - e ) 4 ) (43) 
Using the following change of variable 



we obtain an equation for y 



!=(i+^)v-£ (^) 

(a term a 2 /8 has been neglected in comparison to 1). The solution of ( E^) that satisfies y > (because e > 0) 

depends on the value of the parameter S. One can show that the polynomial P(y) = (1 + -^)y ~\ — § always has a 
positive root. We retain only the smaller positive root y of ( |45| ) because in thermodynamic equilibrium, the system 
always chooses the state with minimal Gibbs free energy. Restoring the usual units, and using ([li]), the nucleation 
fields are found to be: 

„ 00 00 

Hl = yo ^R~\ + ^R 2 

H n+1 =H 1 + (46) 
When the applied field H e lies between H n and H n+ \ 1 the disk contains exactly n vortices and its magnetization is 



calculated using (14). In figure 5, we have plotted the magnetization of a mesoscopic disk with R = 10A\/2 both 
from exact numerical solutions of the Ginzburg-Landau equations and from the expression (|l4|). The agreement 
is very satisfactory. For larger values of the number n of vortices, a discrepancy between the theoretical and the 
numerical expressions appears which results from the interaction between the vortices and the edge currents that we 
have neglected until now. 

The expression ([H]) is also in good agreement with previous experimental and numerical results |l],|7j . A non-linear 
Meissner behaviour still exists before the nucleation of the first vortex as well as between successive jumps. The field 
H\ of nucleation of the first vortex scales as R . The transition between a state with n vortices to a state with 
(n + 1) vortices is of first order since the entrance of a new vortex induces a jump in the magnetization. These jumps 
are of constant height and have a period . If we use the experimental values of for R and A we obtain a value 
for the period of the jumps which is in very good agreement with the experimental value. 

If R is smaller than a threshold value, the system is a fractional fluxoid disk with a second order phase transition. 
If R ~ 1, a vortex can nucleate in the disk and a first order transition occurs. When R increases, the number of 
jumps increases (as R 2 ). These qualitative changes of behaviour with increasing R, which are the important features 
obtained from thepresent model, have been indeed observed in experiments carried out on disks of different sizes. In 
an earlier study |B8[ , we obtained satisfactory values for the nucleation fields but the fractional fluxoid disk, and the 
different regimes obtained by increasing R could not be explained because we neglected subdominant terms that are 
retained here. 

It has been observed experimentally that the period and the height of the jumps cease to be constant when the 
number of vortices increases. These effects are related both to interactions between the vortices and between vortices 
and edge currents. The purpose of the next section is to take into account these interactions and to obtain a better 
estimate for the free energy and the magnetization of a mesoscopic disk. 



VI. WEAKLY INTERACTING VORTICES IN THE VICINITY OF THE DUAL POINT 

So far we have obtained analytical expressions for the free energy and the magnetization of a thin superconducting 
disk at the dual point. When the Ginzburg-Landau parameter has the special value, k — vortices do not interact. 

This fact, discussed in |2j^{|, implies that the bulk free energy does not depend on the location of the vortices. 
However, when k is away from the dual point, the vortices start interacting among themselves; therefore the bulk 
free energy ceases to be a purely topological integer n and the vortex interaction energy must be taken into account. 
Because of this interaction the vortices are no longer necessarily placed at the center of the disk: in an equilibrium 
configuration, the cylindrical symmetry can be broken and the optimal free energy may correspond to geometrical 
patterns such as regular polygons, polygons with a vortex at the center, or even rings of polygons [ ^0p^ , |2q ]. It is the 







competition between the interaction amongst vortices and the interaction between vortices and edge currents that 
determines the shape of the equilibrium configuration. 

Analytical studies were mostly carried out in the limit re — > oo and were based on the London equation ]l7| , p0| , p3t 
for which vortices are point-like and have a hard-core repulsion [fl3|| . We shall study a regime where re is slightly 
different than -^=, i.e. a regime where vortices interact weakly. We shall determine, to the leading order in (re — -^), 
the interaction energy of the vortices. 



A. The interaction energy 

In order to obtain an estimate for the free energy of a system of interacting vortices, we have solved numerically the 
Ginzburg-Landau equations for a cylindrically symmetric infinite system with n vortices located at the center (these 
equations are explicitly written in Appendix A). The free energy per vortex is plotted in figure 6 as a function of re, 
for n = 1,2,3,5 and 10. At the dual point, the free energy per vortex is equal to 1 and is independent of n: all the 
curves pass through this point. When re is different from the interaction between the vortices changes the value 

of the free energy. One can deduce from figure 6 that vortices attract each other for re less than while they repel 

each other when re > -7=. 

From our numerical results we observed that in the vicinity of the dual point, the free energy ^(Kjn) satisfies the 
following scaling behaviour: 

-i-JT( K , n ) = n{ K V2) a{n) (47) 

We note that the relation ( f47| ) is exact at the dual point. For n = 1, 1) is nothing but the self energy Us of a 
vortex. In the vicinity of the dual point we can write: 

i^(/s,l) = l + a(l)(«V^-l) (48) 

The values of the function a(n), as determined from numerical computations, for n ranging from 1 to 30 are given in 
the table | 

We can now derive an approximation for the free energy of a n vortices configuration located at the center of 
the disk and for re close to -^=. Since this configuration is cylindrically symmetric, one can again use the circle T to 
separate the system into two subdomains VL\ and and then estimate separately the two contributions to the total 
free energy. From our numerical scaling result, we deduce a formula for the bulk free energy of a finite system which 
is valid for re close to the dual point. Expanding ( p7[ ) in the vicinity of the dual point, we obtain: 

— T{Q, X ) = n + (kV2 - l)na(n) (49) 
and the boundary contribution, obtained via a variational Ansatz is now given by: 

J_ , s_ ( Av 2 s (R) - B( K )vf(R) if Kn-^ e )|<2re 



where A is still given by the relation (|30|) while B(k) is now given by B(n) — 5/16an 2 . 

The magnetization curve of figure 7 shows both the numerical results and a plot of the magnetization deduced from 
( ^p| ) using (|lj). We notice that the magnetization of a mesoscopic disk is modified when the interactions between 
vortices are taken into account. The period and the amplitude of the jumps are not constant anymore; besides, the 
non-linearity of the curve between two successive jumps is enhanced. These important features of the M — H e curve 
were observed in previous experimental and numerical results Here we have shown that these features are a 

consequence of vortex interactions. 



B. Two-body interaction energy 

The exponent a(n) in the relations (|47]) or ([l9]) allows to describe the interacting potential between vortices. It is 
interesting to compare the result p^ ) with the energy of n vortices obtained by assuming a two-body interaction. In 
this case the energy of the whole system of n vortices can be written as a sum of two terms 
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J_^ = nWs+ <!L_i) Wj( o) (5i) 

Z7T Z 

where Us represents, as noted before, the self-energy of a vortex and Ui the two body interaction potential. Using 
the data of P9fl , we can estimate these two energies to the leading order in — 1). We obtained: 

U s = I + (3i{nV2 - 1) with /3i~0.4 (52) 
Ui(r) = /3 2 (kV2 - 1) min ( 1, cxp(-C(r - -))) (53) 



with 02 — j and C ~ |. From this analysis, and assuming only two-body interaction, we derive an approximate value 
for the free energy of a configuration with n vortices placed at the same point: 



^fr(0) ~ n + ( K V2 - l)n (fa + 



1 nCn — 11 

— F(k, n) = nU s + V ' U^O) ^ n + («V2 - l)n ( fa + fa^^ ) (54) 



If we compare this relation to the previous expression ( p9| ) we find that instead of the sublinear function a(n) we have 
a linear behaviour fa + ^^fa. Hence, the function a(n) takes into account not only two-body interactions among 
vortices but also multiple interactions which are present for values of K around the dual point unlike the large k limit 
where only the two-body contribution remains. 



VII. VORTEX/EDGE INTERACTIONS IN SYSTEM WITHOUT CYLINDRICAL SYMMETRY 

In this section, we calculate the energy at the dual point of a system with only one vortex that is not located at 
the center of the disk. Such a configuration is not in thermodynamic equilibrium and its free energy can be related 
to a surface energy barrier (analogous to the classical Bean-Livingston barrier in the London limit). We first show 
that even when the cylindrical symmetry is broken, the system can still be separated into bulk and edge domains. 



A. Bulk and edge domains. The curve V 



We have seen in section III B that when one or more vortices are located at the center of the disk, there exists a 
circle T on which the current vanishes identically. This circle allowed us to define a bulk and an edge domain and to 
identify the bulk energy with the fluxoid. 

If all the vortices are not placed at the center of the disk (i.e. the configuration is not cylindrically symmetric) there 
is in general no curve of zero current. However the curve T has now the following property: at each point M of T 
the current J*is normal to T. The existence of such a curve is shown by the following argument. Consider a disk with 
only one vortex V situated at a point different from the center of the disk. Take a line segment joining the vortex V 
to the closest point S on the boundary of the disk (see figure 8) . The component of the current density normal to the 
VS segment changes its sign when one goes from V to S. Hence, there exists a point M along this segment where 
the current either vanishes or is parallel to VS. To draw the curve T we start from M in a direction orthogonal to the 
VS segment, and then T is constructed via infinitesimal steps by imposing that at a point M' = M + dM, very close 
to M, the direction of T is orthogonal to the direction of the current at M'. 

Although we lack a general proof, we believe on topological grounds that for vortices at arbitrary positions, there 
always exists a T curve which is everywhere orthogonal to the current (one should note that T does not necessarily 
have only one connected component). In p4j, we shall present a numerical construction of T. In the sequel of this 
work we assume that T exists, that it encircles all the vortices, and consists of one or many simple closed curves. We 
shall call the curve T the separatrix. 

Using r, the domain Q can be decomposed in two regions f2i and such that: 

(i) fiiU^ = fi; 

(ii) fix contains all the vortices ( fii may have multiply connected components); 
(iiii) Q 2 contains the edge of the disk; 

(iv) the separatrix T is the boundary between fl\ and and is everywhere normal to the current density. 
The remarkable property of the separatrix implies that along T one can write: 

(j + A).dl - £ (-^ + A) .dl = j V X .dl (55) 
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since along T, j.dl = 0. Since the separatrix is the boundary of Oi, the property ( |55| ) ensures that the total magnetic 
flux through f^i is quantized. Hence, at the dual point, we can again use the method of Bogomoln'yi and find the 
free energy of f2i to be a purely topological number, just as for an infinite domain, even if the cylindrical symmetry 
is broken. 



B. Free energy of one vortex: the surface energy barrier. 

As before, we estimate the contribution J-^f^) to the total free energy via a variational Ansatz, taking the modulus 
of the order parameter to be constant. To obtain a qualitative result for the surface energy barrier we neglect the 
magnetic energy so that, at the dual point, we have: 

1 -rvn ^ f u.,2,^. T,2 , (1-M 2 ) 2 



JSlo 



2vr 



n 2 

2a (Vo> s 2 ) + (l-^o) 2 ) (56) 



where 



^-\V X -A(R)f (57) 

ITT 



is the superfluid velocity square averaged over the boundary of the disk. As before, we have replaced the integral over 
f^2 by a line integral along the boundary of the sample (i.e. the disk of radius R) multiplied by an effective length 6. 
The function x appearing in (|5^) is the phase of the order parameter, and the vector potential is taken, as before, to 
its value on the boundary of the sample. Optimizing (|56|) with respect to ipo we find that: 

*l = 1 - ^ (58) 

S*fc>- 5 <»> 

for (Vg) < y/2. The phase function \ and the vector potential near the edge of the disk are calculated in Appendix 
B. Using these results, we obtain (for n=l): 

i^(0 2 ) = A (a(l - 4> e f ~ £(1 - <t> e A + f(x, a, </> e - 1)5 (60) 
2ir la \ 4 / 

The function f{x,a 1 (j) e — 1) determines the dependance of the free energy on the position x of the vortex; hence, it 
measures the interaction energy between the edge currents and the vortex as a function of its position. It is given by 

lax 2 , , /. 2 — 1) 



f(x, a, 4> e -l) = Y—ztt* - ^1 ~ « \ _ J ) (61) 

From this expression, we observe that the edge currents tend to confine the vortex inside the system. In figure 9 the 
surface energy as a function of the position x of the vortex is plotted. According to (p|h, only the increasing part 
of the curve is physical. We nevertheless plot the curve defined by ( |6l] ) in the whole range < x < 1 in order to 
emphasize the similarity between our result and the well-known Bean-Livingston surface barrier effect that was first 
derived using the London theory 



VIII. CONCLUSION 



In this work, we have obtained analytical results for the free energy and the magnetization of a mesoscopic super- 
conductor. We have used a known exact solution for the two dimensional Ginzburg-Landau equations in an infinite 
plane, valid at the dual point, to study a finite system with boundaries. With the help of numerical simulations, 
we have carried out a perturbative calculation in the vicinity of the dual point. This approach enabled us to study 
thermodynamically stable states but also metastable states (to obtain a surface energy barrier). This model gives 
theoretical insights into the physical mechanisms involved in the experimental results of |0,H and our analytical results 
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agree quantitatively with experimental measurements. In fact, other related thermodynamic quantities such as the 
surface tension measuring the thermodynamic stability of vortex states can also be computed along this way and 
could generalize to two-dimensional systems previous results obtained in one dimension |36| ]. 

More generally, we believe that a theoretical study in the vicinity of the dual point provides a lot of information 
about the Ginzburg-Landau equations. Although one usually relies on exact results derived from London's equation, 
one should be aware of the fact that these results agree with numerical simulations of Ginzburg-Landau equations 
only when n is large (typically n > 50). We verified that the behaviour we found in the vicinity of the dual point, 
such as the scaling of the free energy, remains valid when k ranges from 0.1 to 10 and this interval of values is indeed 
relevant for many conventional superconductors. 

Our study can be extended in many directions. The scaling results in the vicinity of k = l/v2 were derived from 
numerical simulations: a systematic perturbative expansion around the dual point would put them on a more rigorous 
basis. Secondly, a linear stability analysis of the cylindrically symmetrical solution |35| should allow to understand 
the fragmentation transition between a giant vortex and unit vortices. Since the scparatrix T exists even for vortex 
configurations breaking cylindrical symmetry, our approach can be used to analyze hysteretic behaviour of metastable 
states, and to study polygonal vortex configurations found numerically in mesoscopic superconductors [p|Jlot. 
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X. APPENDIX A: THE GINZBURG-LANDAU EQUATIONS IN A CYLINDRICALLY SYMMETRIC 

SYSTEM 

For a cylindrically symmetric system, we can use ip = f(r)e and A = A(r)eg where n is a non-negative integer 
which represents the number of vortices at the center of the system. We also define the superfluid velocity v s = v s (r)eg, 
where 

«.(r) = - A(r)) (62) 
In this case the Ginzburg-Landau equations are: 

i (~Jrv,)) = 2v,f (64) 



dr \r dr 

It is convenient to define the quantity p(r) — rv s (r). The magnetic field B — B(r)e z is given in terms of p(r) by 

B(r) = --^ (65) 
r dr 

We obtain finally two coupled ordinary differential equations 

f» = -2 K 2 f(l-f)+p*f/r 2 -f/r (66) 
p" = 2pf 2 +p'/r (67) 

with the following boundary conditions at r = a -1 for n^O: 

p{u) = n p(a ) = n — <p e 
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for a disk and 



/(0) = f(a- 1 )=0 ( 

p(0) = n p'ia- 1 ) = -2a(j) e [m) 

for a cylinder. These are the equations we have solved numerically using the relaxation method (38). From the analysis 
of the equations ( |(37j) we deduce the following behaviour in the vicinity of the center of the disk: 

/ ~ r n and p ~ r 2 when r — > 

The free energy (^) is then given in terms of the solution of ( |67j ) by 

J7 /.l/o / B 2 



•2tt 



rdr _ + «*(!-/*) (70) 



XI. APPENDIX B: PHASE AND VECTOR POTENTIAL OF AN OFF CENTERED CONFIGURATION 

WITH ONE VORTEX 

In this appendix we measure the distances in units of R, so the disk has unit radius. Suppose that the vortex is 
located at a distance x from the center of the disk (0 < x < 1). The phase x(p>@) °f the order parameter satisfies 
A% = everywhere on the disk except on the vortex with boundary condition n.V% = 

Using the image method, the phase x(Pi at a point located at a distance p from the center of the disk (with 
< p < 1) is given by @: 

x(p,0)=^{ peMid) _ x - 1 ) (71) 

where Im denotes the imaginary part part of a complex-valued function. Or equivalently: 

/ „ . 1 — x 2 sin 9 
tan X (p,9) = —— 2 (72) 
1 + x cos — .-I 

On the boundary of the disk, p = 1, and one finds that 

d\ 1 — x 2 1 d\ 

~d6 ~ T- YTTgT^^' ^ 

therefore 



= (73) 



g|V x(M )|Wi±!J (74) 

The vector-potential A(R) at the boundary of the sample is a function of the polar angle 9 since the vortex is not 
at the center of the disk. We determine A{R) from the following conditions: 



V.A = 0, * A{R).dl 
Jdn 



and on the boundary A(R).h = 0. The following choice, 

A(R) = e V x (75) 
valid near boundary of the system, satisfies these requirements. 
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FIG. 3. Comparison between B and 1 — \ip\ 2 in a cylinder of radius 10 A\/2 containing one vortex. 
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FIG. 4. Magnetization of a fractional fluxoid disk. Comparison between the experimental measurements [8] (for R = 0.31/im) 
and the theoretical curve taken from the expression (k3) . 
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FIG. 5. Behaviour of the magnetization of a disk with radius 10 A\/2, at the dual point. Dots represent the numerical 
solution and the solid curve the expression (^) together with (|3l|). The only free parameter S has been taken to 8 = 0.76A. 
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FIG. 6. Behaviour of the free energy per vortex F/n = T '/27m as a function of v2« f° r different values of n, the number of 
vortices. At the self-dual point \/2k = 1, the energy Tin) — nT(l) so that the interaction energy between the vortices vanishes 
identically. 
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FIG. 7. Magnetization curve of a disk of radius 10 X\ 
numerical solution and the solid curve the expression (TL 
S = 0.76A. 



/ 2, as a function of the applied field for Ky2 = 0.9. Dots represent the 
j) together with ^). The only free parameter 5 has been taken to 
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FIG. 8. The separation of a system without cylindrical symmetry in two subdomains by a curve V. 
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FIG. 9. Confining energy of a vortex inside a disk due to edge currents. 
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n 


a(n) 


n 


a(n) 


n 


a(n) 


1 


0.417 


11 


0.785 


21 


0.841 


2 


0.544 


12 


0.794 


22 


0.845 


3 


0.613 


13 


0.802 


23 


0.847 


4 


0.658 


14 


0.809 


24 


0.850 





u.oyu 


1 ^ 

10 


n si ^ 

U.olO 


ZO 


n s^ 

U.oOo 


6 


0.715 


16 


0.821 


26 


0.855 


7 


0.734 


17 


0.826 


27 


0.857 


8 


0.750 


18 


0.830 


28 


0.859 


9 


0.764 


19 


0.834 


29 


0.860 


10 


0.775 


20 


0.838 


30 


0.862 



TABLE I. The numerical values of the function a(n) for n ranging from 1 to 30. 
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